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ABSTRACT 

Aims. Accretion rates in low-mass protostars can be highly variable in time. Each accretion burst is accompanied by a temporary 
increase in luminosity, heating up the circumstellar envelope and altering the chemical composition of the gas and dust. This paper 
aims to study such chemical effects and discusses the feasibility of using molecular spectroscopy as a tracer of episodic accretion rates 
and timescales. 

Methods. We simulate a strong accretion burst in a diverse sample of 25 spherical envelope models by increasing the luminosity to 
100 times the observed value. Using a comprehensive gas-grain network, we follow the chemical evolution during the burst and for up 
to 10^ yr after the system returns to quiescence. The resulting abundance profiles are fed into a line radiative transfer code to simulate 
rotational spectra of C**0, HCO^, H'^CO^, and N 2 H+ at a series of time steps. We compare these spectra to observations taken from 
the literature and to previously unpublished data of HCO^ and N 2 H^ 6-5 from the Herschel Space Observatory. 

Results. The bursts are strong enough to evaporate CO throughout the envelope, which in turn enhances the abundance of HCO^ 
and reduces that of N 2 H^. After the burst, it takes 10^-10'^ yr for CO to refreeze and for HCO^ and N 2 H^ to return to normal. The 
H 2 O snowline expands outwards by a factor of ~10 during the burst; afterwards, it contracts again on a timescale of 10^-10^ yr. 
The chemical effects of the burst remain visible in the rotational spectra for as long as 10^ yr after the burst has ended, highlighting 
the importance of considering luminosity variations when analyzing molecular line observations in protostars. The spherical models 
are currently not accurate enough to derive robust timescales from single-dish observations. As follow-up work, we suggest that the 
models be calibrated against spatially resolved observations in order to identify the best tracers to be used for statistically significant 
source samples. 
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1. Introduction 

Stars gain most of their mass while deeply embedded in the 
molecular cloud core out of which they first formed. The ac¬ 
cretion rate of matter from the collapsing core onto the star 
is probably highly variable in time. On cloud scales, turbu¬ 
lence can induce order-of- magnitude var i ations in the accretion 
rates onto individual cores (iPadoan et al.Ll201^ . Within a given 
core, variations of similar magnitude are pos sible from mass¬ 
loading onto an embedde d disk or pseudo-disk (IZhu et alil2009l: 
IVorobvov & Basui l2010t) . The disk becomes unstable at some 
point, resulting in a short-lived accretion burst onto the star. Each 
burst causes a temporary increase in luminosity, heating up the 
circumstellar envelope and altering the chemical composition of 
the gas and dust. This paper analyzes such chemical effects and 
explores how to use molecular lines as a probe of variable accre¬ 
tion rates and timescales. 


Variable or episodic accretion is the leading explanation 
for the wide range of protostella r lumi n osities observed in 


large-scale surveys (iKenvon et al, 

19901: lEvans et all 

200^ 

iKrvukova et alll2012l Fischer et al. 

.l2013llDunham et all 

20i4 


and for EUor a nd EXo r outbursts in more evolved pre-main- 
sequence stars (iHerbigl I1977I; lAudard et al.L l2014t) . Episodic 
accretion can also affect processes such as disk fragmen- 


* Herschel is an ESA space observatory with science instruments pro¬ 
vided by European-led Principal Investigator consortia and with impor¬ 
tant participation from NASA. 


tation ( Stamatelloset^lj^ 201 ll 1201 2l) and lithium depletion 
(iBaraffe & Chabrieii I201C ). It is therefore important to under¬ 
stand the magnitude-frequency distribution and determine how 
often bursts of a certain intensity occur. Based on large-scale 
variability surveys, the strongest bursts - with accretion rates 
of at least 100 tim es the quiescent val ue - happen every 5-50 
kyr per protostar (IScholz et al.L l2013h . The spacings between 
periodic shocks in jets and outflows support quiescent inter¬ 
vals anywhere from 10^ yr down to 10 yr, separating bursts 
with an unknown range of accretion rates (iDevine et al.L 119971: 
iRaga et ^ l2002t lArce et akl l2013l) . High-cadence photomet¬ 
ric surveys show luminosity variations at the 5-50% level on 
times cales of hours to weeks (iBillot et al.L 120121; IStauffer et al.L 
I 2 OI 4 I) . tho ugh only some of that may be related to accretio n 
variability (iMorales-Calderon et al.L 1201 ll: iRebull et al.L l2014l) . 
Despite the incomplete statistics, it appears that stronger accre¬ 
tion bursts occur less frequently than weaker o nes. Such a distri¬ 
bution is supported by n umerical simulations (IZhu et al.L 12009^ 
IVorobvov & BasuLl2010l) . 


Direct observations of luminosity bursts were traditionally 
limited to the relatively evolved Class II sources, where the lack 
of an obscuring dusty envelope makes the bursts visible at op¬ 
tical wavelengths. However, bursts have no w also been detected 
in at least five embedded Class I protostars (lAudard et al.L 201A) 
and e ven in one deeply embedded Class 0 source dSafron et ^ , 
l2015h . The detection of luminosity flares during the earliest 
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phases of star formation lends credence to the notion of episodic 
accretion as a wide-spread phenomenon. 

Accretion bursts have the ability to alter the chemical com¬ 
position of the circumstellar material. For example, the 2008 
burst in EX Lup sparked the production of crystalline silicates 
(lAbraham et al.LlT009l) and b oosted the column de nsities of mid- 
infrared H 2 O and OH lines (iBanzatti et al.L l2012h . In very low- 
luminosity embedded protostars, whose envelopes are too cold 
to produce pure CO 2 ice, the detection of double-peaked 15- 
yum absorption profiles is attrib uted to a higher lumi nosity at 
some unknown po int in the past (iKim et al.ll201 iL 1201 2t see also 
iPoteet et al.ll2013h . 

The chemical effects of embedded accretion bursts primar¬ 
ily result from changes in the temperature of the circumstel¬ 
lar gas and dust. When a protostar enters a burst, the in¬ 
crease i n luminosity can alter th e D/H ratio of water and other 
species (lOwen & Jacaued l2015h . The higher temperatures also 
lead to the evaporation of so me of the icy gra in mantles (iLee 
2007: l^ser & Berginll2012i hereafter Paper I: IVorobvov et al. 
2013h . After the burst ends, the luminosity returns to the quies- 
cent value and the enve lope cools down almost instantaneously 
(I.Iohnstone et al.Ll2013h . However, at typical envelope densities, 
it takes 10^-10^ yr for the gas to freeze back onto the cold dust 
grains. The abundance profiles are out of equilibrium with the 
observed luminosity for all that time. This could explain the 
presence o f spatially extended C' ^O in eight out of 16 low-mass 
protostars (l.le>rgensen et al.Ll2015h . 

The best example to date of episodic accretion chemistry 
in an embedded protostar is the central gap discovered in spa- 
tially resolved observat ions of H'^CO^ J - 4-3 in IRAS 15398 
(iJeirgensen et al.L l2013h . This distribution is consistent with a 
picture where HCO^ and its isotopologs are destroyed by H 2 O 
in the inner envelope when the temperature exceeds 100 K. 
However, the current temperature at the edge of the H'^CO^ gap 
is only ~30 K. The mismatch between the dust temperature and 
the H'^CO^ morphology implies that IRAS 15398 was hotter 
in the past and that the chemistry is still adjusting to the current 
cold environment. Based on the freeze-out timescale of H 2 O, the 
simplest explanation for that hotter past is an accret ion burst that 
happened 100-1000 yr ago (I.Te>rgensen et al.Ll2013h . 

The aim of the current paper is twofold: to explore episodic 
accretion chemistry in a sample of 25 spherical envelope models, 
and to address how certain molecular line ratios can be used as a 
chronometer of when the most recent accretion burst occurred in 
any given source. Section |2] introduces the sample and summa¬ 
rizes the available observations, including new Herschel spectra 
of HCO^ and N 2 H^ 6-5. Section |3] presents the physical and 
chemical models. Section|4]discusses the chemical effects of an 
accretion burst and Sect.|5]shows how to use observed line ratios 
to derive episodic accretion timescales. Lastly, Sect, [^presents 
the main conclusions. 


2. Source sample and observations 

iPaner ll explored the chemical aspects of episodic accretion with 
a set of single-point models at representative densities of 10^, 
10®, and 10^ cm“^. The current paper analyzes the chemistry 
in a sample of observationally determined spherical protostellar 
envelope models, each covering a range of densities and tem¬ 
peratures. We use the predicted abundance profiles to simulate 
molecular line spectra and identify observable diagnostics. 

The envelope models are based on the 29 low-mass em¬ 
bedded protostars targeted in the key prog ram “Water in star- 
forming regions with Herschel” (WISH; Ivan Dishoeck et all 


I 2 OI ll) . Source coordin ates, distances, and othe r basic properties 
are listed in Table 1 of iKristensen et all (l2012h . In order to iden¬ 
tify the best observable diagnostics, we searched the literature 
for observations of molecular lines and ice column densities. The 
most widely available tracers are C'^O J - 2-1, 3-2, and 5-4 
(20, 26, and 15 sources); H'3 cO+ 1-0, 3-2, and 4-3 (18, 17, 
and 16); and N 2 H^ 1-0 (21). In addition, column densities of 
CO 2 , CO, and H 2 O ice are available for 14, 10, and 12 sources. 
Lastly, we present previously unpublished spectra of HCO^ 6-5 
and N 2 H^ 6-5 (16 and 7 sources) obtai ned with the Heterodyne 
Instrument for the Lar-Infrared (HIL L Ide Graauw et all 120101) 
on the Herschel Space Observatory (iPilbratt et al.L 120101) . The 
full data set is summarized in Tables[T]and|2l AnnendixlAloffers 
more details on the observations and the data reduction. 

Lor lack of molecular line observations, we omit NGC1333 
IRAS3, IRAS 12496, R CrA IRS5, and HHIOO IRS from the 
WISH target list. Our final sample the refore contains 25 pro - 
tostars, including IRAS 15398 (Sect. ITl [jprgensen et al.L 1^0131) . 
These sources have bolometric luminosities from 0.8 to 35.7 Lq 
and span the full evol utionary sequ ence from deeply embedded 
Class 0 to late Class I (lLadi[l99l. 


3. Model description 

3.1. Physical framework 


Spherical density and temperature profiles are ayailable from 
iKristensen et al.l (l2012l) for our entire sample of 25 protostars. 
The density follows a power law, «(H 2 ) cc r^P, where the slope 
p is one of three free parameters in a grid of models from the 
one-dimensional continuum radiatiy e trans fer program DUSTY 
(llyezic et al.l Il999t iJdrgensen et al.. l2002l). The o ther two are 
the mass and size of the enyelone. IKristensen et ^ compared the 
model grid output to the obseryed spectral energy distributions 
(SEDs) and submillimeter brightness profiles to find the best-fit 
parameters for each source (see their Table C.l). 

The WISH target list is biased towards bright protostars, 
which are likely to be in a relatiyely actiye phase of accretion. 
Nonetheless, for the purpose of exploring the chemical effects 
of episodic accretion, we treat all sources as if they are cur¬ 
rently in a quiescent phase. Alternatiye assumptions of some or 
all sources currently experiencing a high accretion rate would 
produce a different set of chemical models, but would not alter 
the conclusions in Sect. H] _ 

Within the best-fit models from iKristensen et al.l (l2012h . we 
simulate an accretion burst by changing the stellar luminosity 
and rerunning DUSTY to obtain a new temperature profile. All 
other parameters remain constant. Our goal is to study the chem¬ 
ical effects of a strong burst, such as might happen at most a 
few times during the embedded phase. Specifically, the burst lu¬ 
minosity for each source is set to be 100 times higher than the 
current luminosity; this is the burst intensity inferred for IRAS 


20131) . The simulated bursts last for 100 
20051) . although shorter durations down 


15 398 (iJdrgensen et al 
yr (IVorobypy & Basu . 
to at least 1 yr do not affect our results. After each burst ends, 
we follow the quiescent chemistry for a maximu m of 10® yr to 
coyer the full range of potential timescales from e.g. lScholz et al.l 
(I 2 OI 3 I) . The models are static and ignore any dynamical effects 
of the enyelope collapsing and dis s ipating on ty pical timescales 


of a few 10® yr (lEyans et al.l 1200^1 Visser et all^oTTl) . 

The heating and cooling rates of the dust are fast enough for 
our purposes that the dust te mperature responds ins tantaneously 
to a change in luminosity (iJohnstone et al.L l2013l) . We set the 
gas temperature equal to the dust temperature at all times. As 
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Table 1. Observed integrated intensities (W = f rmbdv' in K km s ') and half-power beam widths (0 in arcsec).“ 


Source 

HCO+ 

6-5 

H‘3C0+ 

1-0 

H‘'CO+ 

3-2 

H"*CO+ 

4-3 

N 2 H+ 

l-O* 

N 2 H+ 6- 

-5" 

References 


W 

6» 

w 

6» 

W 

d 

W 

d 

W 

e 

W 

e 


L1448 MM 

1.37 

43 

2.0 

43 

1.9 

19 

0.92 

21 

12.48 

27 

0.19 

41 

1,2,3,4 

NGC1333 IRAS2A 

1.36 

43 

1.8 

43 

2.1 

19 

2.7 

14 

14.2 

40 

0.41 

41 

1,2,5 

NGC1333 IRAS4A 

2.29 

43 

2.3 

43 

1.4 

19 

0.86 

21 

15.9 

40 

0.78 

41 

1,2,5,6 

NGC1333 IRAS4B 

1.88 

43 

2.1 

43 

0.57 

19 



13.5 

40 

0.64 

41 

1,2,5 

L1527 

1.06 

43 

2.2 

43 

1.1 

19 

0.47 

14 

3.91 

27 

<0.030 

41 

1,2,4 

CedllOIRS4 







0.38 

18 

4.8 

54 



7,8 

BHR71 

0.95 

43 







8.1 

54 



1,8 

IRAS 15398 

0.75 

43 



0.52 

28 



9.4 

54 



8,9 

L483 MM 

1.25 

43 

1.5 

43 

1.7 

19 

1.2 

14 

14.05 

27 

0.29 

41 

1,2,4 

SerSMMl 

4.30 

43 

1.56 

c 

5.8 

19 

3.1 

21 

28 

27 

1.48 

41 

1,5,8,10 

Ser SMM4 

2.08 

43 

1.13 

c 

2.9 

19 

1.2 

21 

49 

27 



1,8,10 

Ser SMM3 

1.74 

43 

<0.1 

c 

1.6 

19 

0.2 

21 

23 

27 



1,8,10 

L723 MM 

0.77 

43 

0.61 

43 

0.94 

19 

0.70 

14 

5.70 

27 



1,2,4 

B335 

0.63 

43 

0.24 

24 

0.46 

28 



6.05 

27 



1,4,11 

L1157 

0.54 

43 

0.89 

29 

0.61 

19 

0.52 

14 

8.38 

27 

0.07 

41 

1,2,4 

L1489 

0.40 

43 

0.96 

43 

0.82 

19 

0.61 

14 

0.26 

40 



2,5 

L1551 IRS5 



1.55 

19 

2.4 

19 

2.2 

14 

15.7 

40 



2,12 

TMRl 



1.1 

43 

0.51 

19 

0.20 

14 

0.37 

40 



2 

TMCIA 



0.44 

19 

<0.27 

19 

<0.13 

14 

7.1 

40 



2,12 

TMCl 



1.26 

19 

<0.27 

19 

<0.12 

14 

5.8 

40 



2,12 

HH46 IRS 







1.1 

18 





7 

Elias 29 

0.28 

43 

0.36 

72 

0.09 

28 







13 

RN091 



0.93 

72 





11 

27 



8,14 


No tes. References: ( 1) t his work; t2'l|.T0rgen sen et al .l2004 (3~l lGregersen et all 1 9^: t4'llEm Drechtinger et al.l 2009l: tSlISan Jose-Garcia et al.l2015 
andiBenz et al.l20 l'a: tOlBlake et alii995l: (7)lyan Kempe n et al.l2 009l: (SllMardonese t alj 1997 1: (9)lGregersen et al.l2000l : flOl lHogerheiide et aT 
119991 tlli lEvans et al.ll2005l: f 12l lOnishi et al.ll2002l: tlSl lBoogert et al.ll2002l: tUl lButner etd! ] | 19951 . Throughout the paper, Class 0 sources 
appear above the horizontal line and Class I sources below it. T^i, is the main-beam temperature, corrected f or beam effic i encies as detailed in the 
references. Typical uncertainties on the integrated intensities are 10-25%. Upper limits are at the 3cr level. See lYildiz et al.l(l201.3h for a compilation 
of C'*0 intensitie s and beam sizes. The N 2 H+ intensities are summed over all hyperfine components (see text). Integrated intensity over a 
20" X 20" region (iHogerheiide et al.Lll999h . 


Table 2. Observed ice column densities {N(X) in 10*^ cm ^). 


Source 

CO 2 

CO 

H 2 O 

References 

L1527 

10 

18 

47 

1 

CedllOIRS4 

12.26 

6.87 


2 

IRAS 15398 

52.16 

8.35 

148.0 

2 

L723 MM 

49.0 



3 

L1489 

16.20 

7.13 

47.0 

2 

L1551 IRS5 

11.86“ 

4.2* 

109.0 

4,5 

TMRl 

4.92“ 


73.8 

4 

TMCIA 

10.98“ 


53.1 

4 

TMCl 

12.71“ 


79.2 

4 

HH46 IRS 

21.58 

7.98 

77.9 

2 

GSS30IRS1 

3.28 

2.53 

15.3 

2 

Elias 29 

6.7 

1.7 

30.4 

6,7 

Oph IRS63 

6.84 

3.49 

20.4 

2 

RN091 

11.66 

7.46 

39.0 

2 


No tes. References: ( l)lAikawa et alj2 0lA (2) Pontoppidan et alJ 2008: 

d; (4 ) IZasowski et al.ll2009): t .5i lTielens et al.(ll991 


(3) 

( 6 ) 


Cook et alJl20Ti1: ty IZasowski et al.ll2009irc 
Boogert et al.ll2000l : (7) iBoogert et al.l 120081 . 


These values 


scaled up by a factor of 1.36 relative to t he published value s from 
IZasowski et al.l (l2009h . as recommended by ICook et ^ 1201 Ih to ac¬ 
count for different CO 2 laboratory band strengths. ® Scaled up by a 
factor of 1.4 to account for different band strengths. 


noted in iPaner 1 the temperature at any point in the envelope 
is proportional to and thus increases by a factor of 3.2 if 
the star becomes 100 times brighter (see also Llohnstone et ^ 


l2013h . As an example, the top panel of Fig.[T]in Sect.|4]shows the 
“quiescent” and “burst” temperature profiles for IRAS 15398. 

3.2. Chemical network 

The basis of our chemical network is the Rate 12 release o f 
the UMIST Database for Astrochemistry (iMcElrov et al.Ll20T^ . 
In addition to standard neutral-neutral and ion-molecule chem¬ 
istry, Rate 12 includes photodissociation and photoionization re¬ 
actions. These are important only at the outer edges of the en¬ 
velopes, which we assume are irradiated by the mean interstellar 
radiation field. The co smic-ray ionizati on rate is set to a standard 
value of 5 X 10 s ' (lDalgarnoll2006l) . X-rays a re not included ; 
they have been detected in one T Tauri FUor (iLiebhart et al.L 
I 2 OI 4 I) . but X-ray luminosities in embedded protostars remain 
highly uncertain and are beyond the scope of this work. 

We expand the network with freeze-out and evaporation 
of all neutral molecules. The freeze-out timescale depends in¬ 
versely on the density and is typically of similar duration (~10^- 
10"' yr) as the quiescent phase between subsequent accretion 
bursts, allowing for an abundance pattern out of equilibrium with 
the stellar luminosity (lLeell200^ iPaner j) . Indeed, it is the long 
freeze-out timescale that drives most of the chemical aspects of 
episodic accretion explored in this work. 

Non-therm al desorp tion via cos mic rays 
( Hasegawa & Herbs!. Il993h or UV radiation (lOberg et al.L 
I2007h is included but only plays a minor role in episodic 
accretion chemistry. Much more important is thermal desorp¬ 
tion. According to spectroscopic observations, protostellar CO 
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evaporates at 25-30 K (iJOrgensen et all l2002l l2004 1201 
lYildiz et'afl l2013h . This is at least 5 K higher than predicted 
from binding energies measured i n the laboratory, espe cially 
th e “pure i ce” value of 855 K from lBisschon et ^ (l2006h used 
in iPaper IL Hence, we now adopt a binding energy of 1307 
K, measured for C O evaporating from amorphous water ice 
(iNoble et all l2012h . Under protostellar conditions, this higher 
binding energy is consistent with the empirical evaporation 
temperature of 25-30 K. Likewise, we increase the binding 
energy of N 2 from 800 K to 1200 K. The binding energy fo r 
atomic O has long been set to 800 K (IWatson & Salpete'Rll972h . 
but new experiments on water ice surfaces offer conclusive 
evidence that this old estimate is too low dM inissa la iMi 
Minissale et al. in prep.). We adopt the current best estimate 
of 1420 K. Th e bindin g energies for other molecules are 
unchang ed from iPaner iL su ch as 5773 K for H 2 O and 2300 K 
forC 02 (iFraser et all 120011 iNoble et akllMl . 

Our network contains two types of grain-surface reactions . 
The f irst is the usual formation of H 2 (iBlack & van Dishoeckl 
Il987h . The second type is simple hydrogenation of C to CH 4 , 
N to NH 3 , and O to H 2 O. This happens one H atom at a time 
at a rate equal to the adsorption rate of H onto the grain surface 
mul tiplied by the relativ e fraction of the reactant molecule in the 
ice (IVisser et alj ^ 20n|. Conversion of CO ice into CO 2 ice, as 
proposed by iKimetalJ (1201111201 2 h . is not included; see Sect. 
l5.1l for details. 

The elemental abundances relative to the total number of H 
atoms are 0.09 fo r He, 1.4 x lOy^ for C , 7.5 x lO^^ for N, and 
3.2 X 10 for O (iMcElrov et al.L l20l^ . At the start of the first 
accretion burst, the ch emical composition is set to typical pre - 
stellar core conditions (iMaret et al.L l 200 (j : IWhittet et ail.L 120091) : 
hydrogen in atomic H (0.005%) and H 2 (~100%); carbon in CO 
(37%), CO ice (35%), and CO 2 ice (28%); remaining oxygen 
in H 2 O ice; and nitrogen in atomic N (55%), N 2 (42%), and 
NH 3 ice (3%). The initial ratio of CO 2 ice to CO ice is 0.8:1, 
the average value observed in three dense molecular clouds be¬ 
lieved to be re presen tative of the earliest stage of star formation 
(IWhittet et aklUPOl . 


3.3. Line radiative transfer 

In order to compare our model to observations, we simu¬ 
late molecular line spectra with the one-dimensional radia - 
tive transfer code RATRAN dHogerheiide & van der TakLl2000li . 
RATRAN solves for the level populations as function of posi¬ 
tion, accounting for both collisional and radiative excitation, and 
then performs ray tracing to synthesize a spectral image. Lastly, 
the images are convolved to the appropriate telescope beam for 
each source and transition (Table [TJ. 

The low-y lines of CO and HCO^ are optically thick for most 
sources, so the optically thin isotopologs C **0 and H'^CO^ 
are used instead. The isotope ratios are set to 69 for '^C/'^C 
and 557 for dWilsonl Il999h. C ollis ion rates ar e take n 

from lFloweddl999t) . lDaniel et al.ld2005l) . and lYang et al.ld 2010 l). 
as compil ed in the Leiden A tomic and Molecular Databasqj 
(LAMDAi l^choier et al.LlTOfBh . 

4. Effects of accretion bursts on abundance profiles 
and rotational spectra 

As described in Sect. 13.11 we simulate an accretion burst in all 25 
sources by increasing the luminosity to 100 times the observed 

* http://home.strw.leidenuniv.nl/~moldata 



Fig. 1. Top panel: radial density profile in IRAS 15398 (black), 
along with the temperatures in the quiescent phase (blue) and 
during the burst (red). Other panels: radial abundance profiles 
of CO, H 2 O, N 2 H+, and HCO^. The red curves are during the 
accretion burst. The black curves are during the quiescent phase, 
at increasing amounts of time after the burst: 10 yr (dotted), 150 
yr (solid), 800 yr (dotted), 2600 yr (short-dashed), 6200 yr (dash- 
dotted), and 24000 yr (long-dashed). The quiescent profiles cul¬ 
minate in the blue curve at 1 x 10 ^ yr after the burst. 


luminosity for a duration of 100 yr. Starting with the pre-stellar 
core composition from Sect. 13.21 we evolve the chemistry dur¬ 
ing the burst at 78-89 radial points per source. After 100 yr, the 
burst ends and the temperature instantaneously returns to the qui¬ 
escent profile. Continuing with the abundances from the end of 
the burst, we evolve the chemistry for another 10 ^ yr and extract 
abundances at a number of intermediate time steps. 

All sources have the same qualitative abundance profiles and 
molecular line spectra, so we choose IRAS 15398 for a quanti- 
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Fig. 2. Synthetic spectra for the quiescent phase in IRAS 15398, scaled as indicated. The different colors correspond to the same 
post-burst time steps as in Fig. [1] 10 yr (black), 150 yr (red), 800 yr (gray), 2600 yr (purple), 6200 yr (green), 24000 yr (orange), 
and 1 X 10^ yr (blue). The arrows indicate how the intensity of each line changes as function of time; the colored flags mark the 
peak intensities at the listed times. The 1-0 panel shows the isolated Fi = 0-1 hyperflne component, whereas the 6-5 panel 
shows the sum over all unresolved components. 


tative discussion. The results for the rest of the sample are pre¬ 
sented in Sect.|5] 


4. 1. Abundance profiles 

IRAS 15398 currentl y has a luminosity of 1.6 Lq and an enve¬ 
lope mass of 0.5 Mq (iKristensen et al.L[2012h . The top panel of 
Fig.IUplots the temperature profiles for the quiescent phase (Tq) 
and the burst (Tb) in IRAS 15398, along with the density profile 
(«). The quiescent temperature decreases from 250 K at the inner 
edge of the envelope to 10 K at the outer edge. During the burst, 
the entire curve goes up by a factor of 3.2 (Sect.[3Tj. 

The other panels of Fig.[T]show the radial abundance profiles 
of CO, H 2 O, N 2 H^, and HCO^ during the burst (red curves) and 
at a series of time steps after the burst (black and blue) . The basic 
results for CO are the same a s reported elsewhere (lLeell200^ 
iPaner IL IVorobvov et al.ll20T^ . The temperature exceeds 30 K 
everywhere during the burst and more than 99.9% of CO is in 


the gas phase at all radii0 After the burst, CO freezes out beyond 
250 AU, where the quiescent temperature lies below 25 K. Since 
the freeze-out rate is proportional to the gas density, depletion 
is fastest right around the 25 K radius. For the specific density 
profile in the envelope model of IRAS 15398, it takes 700 yr for 
the CO abundance to drop by a factor of 10 and 1500 yr for a 
factor of 100. 

The column density of CO ice integrated through the enve¬ 
lope is reduced almost to zero during the burst. Once the en¬ 
velope cools down, the ice column starts to reform. It reaches 
20% of the original amount after 100 yr and 90% after 4000 yr. 
CO 2 has a higher binding energy than CO and evaporates around 
45 K. This allows 15% of the integrated CO 2 ice column to sur¬ 
vive during the simulated burst in IRAS 15398. It takes 600 yr 
in the quiescent phase to rebuild 90% of the original CO 2 ice. 

H 2 O evaporates out to 100 AU during the burst. Like CO, 
the excess H 2 O refreezes in the quiescent phase, but it does so 


^ I Led (l2007h simulated a lower-luminosity protostar and retained 
some CO ice in the outer envelope during the burst. 
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about ten times faster because of the higher densities in the inner 
envelope. All H 2 O that evaporated during the burst disappears 
within a few 100 yr in the envelope model of IRAS 15398. On 
timescales of a few 1000 yr, the freeze-out of CO removes a 
substantial amount of oxygen from the gas phase. This leads to a 
drop in the H 2 O profile near the 25 K radius at 250 AU. The H 2 O 
abundance increases again towards the outer envelope because 
of photodesorption. The integrated column density of H 2 O ice 
drops by 65% during the burst; afterwards, it takes only ~10 yr 
to rebuild the column to 90% of the equilibrium value. 

The chemistry of N 2 H^ is strongly tied to that of CO 
and H 2 O, because both act as important destruction channels. 
The depletion of CO and H 2 O in the quiescent phase thus al¬ 
lows N 2 H^ to increase in abundance at almost all radii for the 
first few 1000 yr after the burst. This type of anti-correlation 
between CO and N 2 H^ is well known from observations of 
molecular clouds, embedded protostars, and circumste llar disks 
(iBergin et al.L 12002 : J^gensenL 12004 lOi et al.L 1201 3h and was 
also noted bv iLeel ( 2007h . On timescales of more than a few 1000 
yr, freeze-out of N 2 near 400 AU causes a decrease in the N 2 H^ 
abundance. 

Lastly, the abundance of HCO^ undergoes a seesaw mo¬ 
tion around the 25 K radius. CO is the primary parent species 
of HCO^ and H 2 O is one of the dominant destroyers. Freeze- 
out of H 2 O allows more HCO^ to be formed in the inner en¬ 
velope, while freeze-out of CO reduces the amount of HCO^ 
in the outer envelope on longer timescales. The correlation be¬ 
tween the protostellar abundances of CO and HCO^ is ano ther 
well-established observational result (Ijprgensen et alil2004ll . 

The chemical evolution in the other 24 envelope models is 
qualitatively the same as in IRAS 15398, but the timescales on 
which the aforementioned changes occur differ depending on the 
density profiles. IRAS 15398 has relatively high densities in the 
H 2 O and CO freeze-out zones compared to the rest of the sam¬ 
ple. As a consequence, freeze-out in the other sources is typically 
slower by up to two orders of magnitude and the chemical effects 
of the accretion burst take longer to subside. 

In summary, the chemical composition of a protostellar en¬ 
velope changes substantially during an accretion burst. These 
changes last for long after the system returns to quiescence: 


• enhanced CO in the outer envelope for a few 10^-10^ yr; 

• reduced CO ice and CO 2 ice in the outer envelope for a few 
10 ^-10^ yr; 

• enhanced H 2 O in the inner envelope for a few 10^-10^ yr; 

• reduced N 2 H^ in the outer envelope for a few 10^-10"^ yr; 

• seesaw pattern for HCO^: enhanced in the inner envelope, 
reduced in the outer envelope for a few lO^-lO'^ yr. 


The exact timescales in a given source depend on the density 
profile. The next section explores the consequences of the abun¬ 
dance changes on the rotational emission lines. 


4.2. Rotational spectra and ice column densities 

The abundance changes caused by an accretion burst persist for 
long after the burst has ended. To what extent can the chemical 
effects of a burst still be observed once the protostar has returned 
to quiescence? Figure |2] shows a set of nine single-dish spectra 
corresponding to the abundance profiles at the quiescent time 
steps from Fig.[T] All spectra are first convolved with the appro¬ 
priate telescope beam (Table [TJ and then continuum-subtracted. 
The arrow in each panel indicates schematically how the peak 
intensity changes with time: continuous decrease, continuous in¬ 
crease, or an increase followed by a decrease. 



Fig. 3. Simulated integrated line intensities and ice column den¬ 
sities in IRAS 15398 as function of time during the quiescent 
phase. The N 2 H^ J = 1-0 and 6-5 intensities are summed over 
all HF components. Where available, the gray bars show the ob¬ 
served intensities and columns from Tables [1] and with 20% 
uncertainties. 


For a more quantitative view of the temporal behavior, the 
spectra from Fig.|2]are integrated and plotted as function of time 
in Fig. [3] Also shown here are the simulated column densities 
of CO 2 , CO, and H 2 O ice. Where available. Fig. [3] includes the 
observed intensities and column densities from Sect. |2] as gray 
bars with 20% error margins. We compare the model results to 
the observations in Sect. EU the remainder of the current section 
discusses the model trends on their own. 

During the quiescent phase, the simulated line intensities 
change predictably based on the abundance variations. The 
freeze-out of CO leads to a monotonic decrease in the optically 
thin C'^O lines. H'^CO^ 1-0 and 3-2 are dominated by the cold 
outer envelope and their intensities also decrease monotonically 
as CO and HCO^ are depleted. H'^CO^ 4-3 and HCO^ 6-5 orig¬ 
inate closer to the star and pick up on the seesaw abundance 
pattern from Fig. [1] the line intensities show an initial rise be¬ 
fore depletion in the outer envelope sets in and the intensities go 
down. 

Because of the abundance anticorrelation between CO and 
N 2 H^, the line intensities of C'^O and N 2 H^ evolve largely in 
opposite directions. As CO and H 2 O freeze out after the burst, 
N 2 H^ becomes more abundant and its 1-0 and 6-5 lines gain in- 
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tensity. The 1-0 intensity decreases again on timescales of more 
than ~24 000 yr after the burst. This is due to the depletion of 
atomic N and the resulting loss of N 2 H^ in the cold outer enve¬ 
lope. The 6-5 line is not very sensitive to the N 2 H^ abundance 
around 1000 AU; its intensity levels off at late times, but does 
not turn over. 

The total observable column densities of CO 2 and CO ice 
increase by factors of 3 and 20 on timescales of 10^-10"^ yr after 
the burst. The column of H 2 O ice is only marginally affected by 
changes in luminosity and remains constant at all times. 

5. Line ratios as a diagnostic of episodic accretion 
timescaies 

One of the challenges in episodic accretion is to understand the 
magnitude-frequency distribution; how often do burst of a cer¬ 
tain magnitude or intensity occur? The chemical effects explored 
in this work can help constrain the frequency of very strong 
bursts, where the accretion rate increases by at least a factor 
of 100. In the quiescent phase following such a burst, several 
observables vary monotonically with time. In principle, any of 
these can be used as a chronometer of when the most recent burst 
occurred. For a statistical sample of protostars, it then becomes 
possible to determine the average burst frequency. 

With the current set of spherical models, inferring when the 
most recent burst occurred in any given source will be a rather 
crude method with uncertainties at the order-of-magnitude level. 
The following discussion should therefore be seen as an attempt 
to explore certain methodologies, rather than an attempt to de¬ 
rive firm numbers. For an alternativ e approach based on sp atially 
resolved observations of C^*0, see lJOreensen et all (l2015h . 

5.1. Models vs. observations 

Figure [3] highlights both the successes and the shortcomings of 
our model for IRAS 15398. Six of the observables are repro¬ 
duced at some point during the quiescent phase, but the times 
at which the matches occur range over more than two orders of 
magnitude: from ~30 yr for CO ice to ~10"^ yr for C'^O 2-1. 
Furthermore, H 2 O ice is overproduced at all times, and HCO"^ 
6-5 and N 2 H^ 1-0 are underproduced. 

The overproduction of H 2 O ice may be due to a lack of 
grain-surface reactions in our chemical network, which tend 
to drain oxygen out o f H 2 O and convert it to other species 
dSchmalzl et al.L I20l4) . N 2 H^ 1-0 probably suffers from our 
choice of a spherically symmetric envelope model. In spatially 
resolved observations, the morphology of N 2 H^ 1-0 often de¬ 
viates significantly from spherical symmetry (I.l0rgenseni l2004t 
iTobin et all l20 1 ih . This can result in enhanced emission in one 
direction, e.g. from interactions between the outflow and the 
dense envelope, which our spherical models do not reproduce. 
The underproduction of HCO^ 6-5 appears to be an excitation 
effect; this line is enhanced along the outflow cavity walls due 
to direct ul traviolet heating, similar to what is see n in *^CO and 
6-5 (ISnaans et al.LI199.5HYddiz et aTll2ni2h . 

iKim et aP (l20 1 iL |20 1 2h argued that some fraction of CO ice 
needs to be converted into CO 2 ice during each quiescent phase 
in between subsequent bursts. However, that conclusion may be 
an artifact of the low initial CO 2 abundance in their model. Our 
model starts with the avera ge CO and CO 2 ice a bundances ob¬ 
served towards dense cores (IWhittet et al.Ll20(M . As a result, it 
does not require any conversion of CO into CO 2 ice to reproduce 
the observed column densities in IRAS 15398. 


Table 3. Candidate observable ratios to be used as chronometers 
of episodic accretion timescales. 


C‘'*0 2-l /N 2 H+ 1-0 

C‘*0 3-2/N2H+ 1-0 

H‘3CO+ 1-0 / N 2 H+ 1-0 

H'3CO+ 3-2 / N 2 H+ 1-0 

C‘®0 5-4 / N 2 H+ 6-5 

C‘**0 2-1 / CO 2 ice 

C“*0 2-1 / CO ice 

C‘®0 3-2/C02 ice 

C'*0 3-2 / CO ice 

H‘3C0+ 1-0 / CO 2 ice 

H‘3C0+ 1-0 / CO ice 

H‘3CO+ 3-2 / CO 2 ice 

H‘3CO+ 3-2 / CO ice 


The comparison between model predictions and observa¬ 
tions is expanded to all 25 protostars in Fig. |4] The trends iden¬ 
tified in Figs.|2]and[3]for IRAS 15398 generally hold for the full 
sample: all C'^O lines become weaker as function of time and 
usually match the observed intensities at some point. H'^CO^ 

1- 0 and 3-2 also decrease monotonically and are in reasonable 
agreement with the observations. H'^CO^ 4-3 and HCO^ 6-5 
reflect the seesaw abundance pattern of HCO^ by first increasing 
and then decreasing in strength. Both lines are too weak in the 
models by a factor of 2-10. The N 2 H^ 1-0 line gains intensity 
during the quiescent phase in all sources and turns over at late 
times in about a third of the sample. N 2 H^ 6-5 initially increases 
and then levels off. Both N 2 H^ lines tend to be underproduced 
relative to the observed intensities. The three ice column densi¬ 
ties always increase with time after the burst. CO 2 ice and CO 
ice usually match the observations, but H 2 O ice is consistently 
too abundant in the models. 

In order to construct a consistent set of models, we treated 
all sources as if they are currently in the quiescent phase (Sect. 
HU- This is undoubtedly a false assumption - especially for the 
brighter sources in our target list - and adds another source of 
discrepancies between models and observations. Future stud¬ 
ies of individual sources or large samples will benefit from 
independent constraints on the accretion rates, e.g. from hy¬ 
drogen Bry (IConnellev & Greeii^. l2010l) . If the stellar mass 
is known from me asurements of a rotationally supported disk 
(iTobin et al.L 12012h . the accretion rate can als o be computed 
from the relationship between Lboi, M,, and M (lAdams & ShuL 

EMI). 

5.2. Line ratios 

The absolute line intensities in Fig. [3] are sensitive to uncertain¬ 
ties in the source distance and in various model parameters like 
the envelope mass and the temperature profile. A better approach 
is to take the ratio between two rotational lines whose intensi¬ 
ties change in opposite directions after the burst. This effectively 
acts as a self-normalization against distance and source model, 
in particular for line pairs at similar excitation levels. 

Table a lists thirteen observable ratios as candidate 
chronometers, chosen in part for their sensitivity to the burst, and 
in part for being easily observable with modern sub-millimeter 
facilities. The first four ratios are the low-T line pairs of C'^O 

2- 1, C'^O 3-2, H'3cO+ 1-0, and H'3 cO+ 3-2 overN 2 H+ 1-0. 

The flfth one is the high-7 pair of C'*0 5^ over N 2 H+ 6-5. The 
remaining eight are the ratios between the intensities of C'^O 
2-1, 3-2, H'3C0+ 1-0, and 3-2 on the one hand 

and the column densities of CO 2 or CO ice on the other hand. 

As an example, the left panel of Fig. |5] shows how the ra¬ 
tio between the integrated intensity of C**0 2-1 and the column 
density of CO ice evolves as function of time after the burst for 
the full source sample. The ratio plots for the other 12 pairs of 
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Fig. 4. Comparison between model predictions and observations. This plot is similar to Fig.[3l except that the model values are now 
divided by the observed quantities. A ratio of unity means the model is in exact agreement with observations. If the ratio is larger 
than unity, the model overpredicts the observation, and vice versa. The horizontal scale for each source runs logarithmically from 
10 to 10^ yr after the burst. The observable quantity in the bottom three panels is the column density of ice along a pencil beam to 
the central star. In the other panels, it is the integrated intensity of various rotational lines. The N 2 H^ J — 1-0 and 6-5 intensities 
are summed over all HF components. The sources are ordered from left to right by increasing luminosity. 


observables in Table [3] look qualitatively the same. Each curve 
represents one source model. The ratios decreases with time for 
all sources, but the curves are scattered across more than four or¬ 
ders of magnitude and are of little diagnostic value in this form. 

The scatter in Fig. |5] is due to the range of envelope pa¬ 
rameters encountered across the sample. CO freezes out at the 
same temperature in each source (< 25 K), but because of 
all the different envelope parameters, it does not do so at the 


same density. In high-density envelopes, the freeze-out timescale 
is shorter and the CO gas/ice abundance ratio decreases more 
rapidly than in lower-density sources. Likewise, the C'^0/N2H^ 
and H'^C 0 ^/N 2 H^ line intensity ratios change at different rates 
for different sources. 

The scatter in the left panel of Fig. |5]can be reduced by ex¬ 
pressing the time axis for each source in units of its characteristic 
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Fig. 5. Left: ratio between the integrated intensity of C'^O 2-1 and the column density of CO ice as function of time during the 
quiescent phase. Each gray curve corresponds to one source. Right: same, but with the time axis normalized by the characteristic 
freeze-out timescale for each source. 


freeze-out timescale (iPaner J. erratum): 


Tfr = 1 X lO'^yr 


10 K 10®cm-3 

Efr nfi- 


( 1 ) 


The characteristic freeze-out density is defined as the density 
at which the temperature reaches Tfr in the quiescent phase. We 
tested Tfr from 10 to 40 K and found the greatest reduction in 
scatter for 15 K, where the bulk of the low-7 C**0 emission 
originates. At 15 K, ranges from 3.5 x 10^ cm“^ in Elias 29 
to 5.5 X 10® cm“^ in Ser SMM4. The corresponding freeze-out 
timescales vary from 2.3 Myr to 1500 yr. 

The time-normalized C'*0 2-1/CO ice ratios appear as the 
gray curves in the right panel of Eig.|5] The sources are still scat¬ 
tered across an order of magnitude, but all follow the same trend 
and an average relationship between normalized post-burst time 
and gas/ice observable ratio can easily be drawn. Our source 
sample is probably diverse enough (albeit not unbiased) that we 
can use this average relationship also for other samples where 
detailed envelope models are not readily available. 


5.3. Future prospects 

Using all the observable ratios from Table |3l we can estimate 
when the most recent accretion burst occurred in each of our 
sources. However, the results are currently too scattered and the 
uncertainties too la rge to produce reliable numbers worth tabu¬ 
lating. Eor example, Jprgensen et al.l (1201 3h concluded from spa¬ 
tially resolved H'^CO^ 4-3 that IRAS 15398 experienced a burst 
between 100 and 1000 yr ago. We find the same timescale from 
the C*®0 gas/COa ice ratios, but C**0 gas/CO ice suggests a 
longer timescale of ~6000 yr, whereas H'^CO^ 3-2/N2H^ 1-0 
yields a value as long as 5 x 10"^ yr. Eurthermore, the observed 
C'^ 0 /N 2 H^ line ratios do not match any timescale from the 
model. The situation is no better for most of the other sources. 
This lack of consistency between different tracers was already 
evident in Eigs.[3]and|4l where some curves match the observa¬ 
tions at very different post-burst times (Sect. lSTI) . 

The most likely reason for these discrepancies is the as¬ 
sumption of spherical symmetry for the envelope models (Sect. 
im i. However, even if the spherical models do not produce re¬ 
liable timescales, they are still useful as a basis for follow¬ 


up studies. Spatially resolved observations offer many advan- 
tages over single-dish data (as explored for e.g. C'^O 2-1 by 
iJdrgensen et alJ 1201 5h and can be carried out for statistically 
significant source samples with current facilities. We will ex¬ 
pand the models beyond spherical symmetry and calibrate them 
against existing interferometric data to identify the most reliable 
tracers of episodic accretion timescales. These tracers can then 
be used to derive timescales for many tens or hundreds of proto¬ 
stars and compared against results from other methods. Working 
with such large samples also reduces the problem that any given 
source may be an outlier. 


6. Conclusions 

This paper simulates the chemical effects of episodic accretion 
in a diverse sample of 25 protostellar envelope models. The in¬ 
tensity of the accretion bursts is set to 100 times the quiescent 
luminosity, enough to evaporate CO ice all the way to the outer 
edge of the envelope in each source (Sect. IrTI) . The H 2 O snow¬ 
line expands in size by a factor of ~10 during the burst. The 
evaporation of CO and H 2 O leads to the destruction of N 2 H^ 
throughout the envelope. HCO^ is enhanced outside of the H 2 O 
snowline and destroyed inside of it. 

When the burst ends, the temperatures quickly return to nor¬ 
mal. The species that evaporated during the burst start to refreeze 
onto the cold dust, but this is a relatively slow process. H 2 O re¬ 
mains enhanced for 10^-10^ yr and CO for an order of magni¬ 
tude longer. The abundances of N 2 H^ and HCO^ are reset on 
similarly long timescales. 

The abundance changes resulting from an accretion burst 
also have a substantial effect on various molecular lines, as quan¬ 
tified by radiative transfer simulations at a series of time steps 
(Sect.EJli. C'*^0 2-1, 3-2, and 5-4 and H'^CO^ 1-0 and 3-2 
are enhanced during the burst and decay monotonically after the 
burst ends. H'^CO^ 4-3, HCO+ 6-5, and N 2 H+ 1-0 first increase 
after the burst and then decrease. The column densities of CO 
and CO 2 ice increase monotonically after the burst. 

The simulated line intensities and ice column densities are 
compared to single-dish observations gathered from the litera¬ 
ture and to previously unpublished spectra of HCO^ and N 2 H^ 
6-5 from the Herschel Space Observatory. Several line ratios 
are identified that can be used, in principle, as a chronometer of 
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when the most recent strong burst occurred in each source (Sect. 
|5]l. In practice, however, the spherical source models are not 
accurate enough to derive reliable timescales from single-dish 
data. Interferometric observations are therefore recommended 
for follow-up work. In order to tackle a statistically significant 
source sample, the models need to be calibrated against spatially 
resolved observations to identify robust tracers of episodic ac¬ 
cretion timescales. 
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Appendix A: Details of observations 

The spectroscopic observations used in this work are summa¬ 
rized in Tables [1] and |3 The following subsections offer more 
details on how the data were obtained and reduced. The HCO^ 
and N 2 H^ 6-5 spectra are published here for the first time; all 
other observations were gathered from the literature. Typical cal¬ 
ibration uncertainties are between 10% and 25%. 
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A.1. C‘^0 

lYildiz et^ (l2013h compiled all available C'^O data from J - 
2-1 up to 10-9, with one omission; an integrated intensity of 
ll.OKkms-' in a 22.9" beam for the 2-1 line in Ser SMM4 
dHogerheiide et al.lll999h . The 6-5 and higher lines are detected 
in only a handful of sources and are therefore of limited use as an 
observational diagnostic. The 2-1, 3-2, and 5^ lines have de¬ 
tection rates of 77%, 100%, and 58% and are well suited for our 
purpose. Furthermore, with upper-level energies of 16, 32, and 
79 K, these three lines together are a g ood probe of the a bun- 
dance jump at the CO sublimation front (lYildiz et al.L 120131) . 

A.2. HCO+ 

The low-7 lines of HCO^ are often optically thick, so we 
only consider the 6-5 line at 535 GHz (Fu/^ = 90 K). This 
transition was observed with Herschel-HIFI as part of WISH 
and the related open-time program “Water in low-mass proto¬ 
stars; the William Herschel line legacy” (WILL; proposal code 
OT2_evandish_4; Mottram et al. in prep . ). Fiv e of the spectra 
were p resented bv ISan .lose-Garcfa et al.l (120 1 5h and iBenz et alJ 
(l20I5h . The other eleven HCO^ lines from this dataset are as yet 
unpublished. 

Table lA.ll lists the observing dates and identification num¬ 
bers (ObsIDs), including several duplicate or triplicate ob¬ 
servations taken in spectral setups of other lines. Five of 
these come from the open-time program “Searching for the 
onset of energetic particle irradiation in Class 0 protostars” 
(OT2_cceccare_4). 

All observations were taken in double beam switch mode. 
We use the spectra recorded with the wide-band spectrometer 
(WBS; 1.1 MHz reso lution) and perform basic standard process¬ 
ing in HIPE vlO.0.0 (IOtlll201(ih . followed by further analysis in 
CLASifl The H and V polarization spectra were averaged after 
individual inspection and the intensities were converted to main- 
beam temperature scale through a main-beam efficiency of 0.72 
(iRoelfsema et al.L l2012h . Where available, the multiple epochs 
of data all match to within the calibration uncertainty of 10% 
and are averaged to reduce the noise. This results in final rms 
noise levels of 3-7 mK in 0.28 km s"' bins. The final data re¬ 
duction step was to subtract a linear baseline. 

In all but two sources, the HCO^ 6-5 spectrum shows two 
emission features; a broad component with a full width at half 
maximum (FWHM) of 4-14 km s ' and a narrow component of 
1-2 km s *. These components are also detected in various CO 
lines and correspon d to the out fl ow (b r oad) and the quiescent 
envel ope (narrow; lYildiz et"^ 120131 ; ISan .lose-Garcfa et al.L 
l2015h . The two exceptions in our sample are L723 and B335, 
which only have the narrow component. 

Since we are interested exclusively in the quiescent enve¬ 
lope, we fit two Gaussian profiles to each HCO^ spectrum (ex¬ 
cept L723 and B335) and subtract the outflow emission. Figure 
I A. II sho ws the resulting spec tra, corrected for the source veloci¬ 
ties from lYildiz et al.l|2012h . Table[T]lists the integrated intensi¬ 
ties for the full sample. 

A.3. 

Complementing the HCO^ 6-5 line are the 1-0, 3-2, and 4-3 
lines of optically thin H'^CO"'". Their upper-level energies are 4, 


^ http://www.iram.fr/IRAMFR/GILDAS 


Table A.l. Herschel observations of HCO^ and N 2 H^. 


Source 


ObsID 

Date 

Ref." 


HCO+ 

6-5 at 535.062 GHz 


L1448 MM 


1342203186 

2010-08-18 

1 

NGC1333 IRAS2A 

1342192206 

2010-03-14 

1 



1342202024 

2010-08-01 

1 



1342248912 

2012-07-31 

2 

NGC1333 IRAS4A 

1342192207 

2010-03-15 

1 



1342202023 

2010-08-01 

1 



1342248914 

2012-07-31 

2 

NGC1333 IRAS4B 

1342192208 

2010-03-15 

1 



1342202022 

2010-08-01 

1 



1342202033 

2010-08-02 

1 

L1527 


1342203188 

2010-08-18 

1 



1342250193 

2012-08-24 

2 

BHR71 


1342200755 

2010-07-07 

1 

IRAS 15398 


1342266008 

2013-03-05 

3 

L483 MM 


1342207582 

2010-10-28 

1 

Ser SMMl 


1342194463 

2010-04-11 

1 



1342207581 

2010-10-28 

1 



1342251637 

2012-09-29 

2 

Ser SMM4 


1342194464 

2010-04-11 

1 

Ser SMM3 


1342207580 

2010-10-28 

1 

L723 MM 


1342219172 

2011-04-21 

1 

B335 


1342219182 

2011-04-21 

1 

L1157 


1342199077 

2010-06-23 

1 



1342246462 

2012-05-17 

2 

L1489 


1342203187 

2010-08-18 

1 

Elias 29 


1342266143 

2013-03-07 

1 


N 2 H+ 

6-5 at 558.967 GHz 


L1448 MM 


1342238644 

2012-02-03 

4 

NGC1333 IRAS2A 

1342225935 

2011-08-09 

4 



1342248913 

2012-07-31 

2 

NGC1333 IRAS4A 

1342238646 

2012-02-03 

4 



1342248915 

2012-07-31 

2 

NGC1333 IRAS4B 

1342238645 

2012-02-03 

4 

L1527 


1342250194 

2012-08-24 

2 

L483 MM 


1342244045 

2012-04-10 

4 

Ser SMMl 


1342230370 

2011-10-08 

4 



1342251636 

2012-09-29 

2 

L1157 


1342246461 

2012-05-17 

2 


Notes. References; (1) WISH; [San .lose-Garcfa et ^1201.51 . iBenz et all 
120151 ; (2) OT2_cceccare_4; (3) OT2_evandish_4; Mottram et al. in prep.; 
(4) OTLphilybla.1. 


25, and 42 K, again offering a good probe of chemical changes 
induced by the evaporation or freeze-out of CO. 

The integrated intensities compiled in Table [T| originate 
mostly fro m the ground-based s urveys of iHogerheiide et al.l 
(1199911 and LlOrgensen et all (l2004ll . with other references as in¬ 
dicated. Some sources have been targeted with multiple tele¬ 
scopes, in which case we choose the observation with the small¬ 
est single-dish beam. 

A.4. A/2H+ 

The final two molecular lines are the 7 = 1-0 and 6-5 transi¬ 
tions of N 2 H^ (Eii/ k = 4 and 94 K), bo th of which are prone to 
hyperflne splitting (ICaselli et al.L Il995h . The quantum numbers 
for the HF levels are F\ and F. The 1-0 transition breaks up 
into one isolated fine {F\ - 0-1) and two clusters of three lines 
{F\ - 1-1 and 2-1), which are readily resolved with heterodyne 
spectrographs. The 6-5 transition is dominated by nine HF com- 
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Fig.A.l. Spectra of HCO^ 6-5, corrected for source velocity 
and scaled as indicated. If a broad emission component from the 
outflow was present (in all sources but L723 and B335), it was 
fitted with a Gaussian profile and subtracted from the spectrum 
to leave only the narrow emission from the quiescent envelope. 



Fig. A.2. Spectra of N 2 H^ 6-5, corrected for source velocity and 
scaled as indicated. 


contribution from the outflow. Table [T] lists the integrated inten¬ 
sities for the full sample, including the upper limit for L1527. 


A. 5. Ices 

The ice column densities in Tablel^are based on published mid- 
infrared observations from the ground and from space. The con¬ 
version from spectrum to column density requires an absorption 
band strength measured in the laborat ory for the appropriate pure 
or mixed ices (iGerakines et al.Lfl995h . Most papers on protostel- 
lar ices use the same band strengths and we adopt the column 
de nsities as publis h ed. Th e two exceptions are the CO?, columns 
of Zasowski et al.l (l2009l) and the CO column of iTielens et al.l 
(Il991h . which we scale by a common factor of 1.4 to be consis¬ 
tent with other studies. 


ponents within 0.1 km s ' and therefore appears as a single line 
in astronomical observations. 

Most of the J - 1 -0 in tensities in Table [Ha r e taken from 
iJdrgensenet akl (l2004l) and lEmprechtinger etaH (l2009h . who 
provided integrated intensities summed over all HF co mponents. 
Additional data are taken from iMardones et ^ (Il997h . who only 
showed the spectra of the isolated component and did not tabu¬ 
late any intensities. We measure the integrated intensities from 
their Figures 1 and 2 and multiply by a factor of 9 (based on the 
Finstein A coefficients and level degeneracies) to approximate 
the intensity for the full 1-0 band. 

The 6-5 lines were observed primarily as part of the 
Herschel open-time program “The chemistry of nitrogen in dark 
clouds” (OTl_philybla_l), with two additional sources (F1527 
and FI 157) covered by OT2_cceccare_4. Both programs targeted 
N 2 H+ 6-5 in NGC1333 IRAS2A and 4A and Ser SMMl. None 
of these spectra have been published, so we follow the same pro¬ 
cedure as outlined for HCO^ 6-5 in Sect. lA.2l Table lA.ll lists the 
observing dates and ID codes and Fig. IA.2I shows the reduced 
spectra. The rms noise is 10-14 mK in 0.27 km s“' bins. N 2 H^ 
6-5 is not detected in F1527. All other sources show a single 
narrow emission feature with an FWHM of 1-2 km s ' and no 
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